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Abstract 

We study the influence of repulsive interactions on an ensemble of coupled excitable 
rotators. We find that a moderate fraction of repulsive interactions can trigger global firing 
of the ensemble. The regime of global firing, however, is suppressed in sufficiently large 
systems if the network of repulsive interactions is fully random, due to self-averaging in its 
degree distribution. We thus introduce a model of partially random networks with a broad 
degree distribution, where self-averaging due to size growth is absent. In this case, the regime 
of global firing persists for large sizes. Our results extend previous work on the constructive 
effects of diversity in the collective dynamics of complex systems. 

1 Introduction 

Noise and disorder can have a variety of effects on the collective dynamics of complex sys- 
tems such as biological populations, chemical reactions, oscillator ensembles, among others. 
Somehow paradoxically, such effects often play a constructive role in inducing coherent be- 
haviour in the system [Hill [3]. An already classical example is the phenomenon of stochastic 
resonance, where noise of appropriate intensity enhances the response of a system to an os- 
cillatory external forcing [5. A similar effect in an extended system can be caused by the 
presence of some degree of heterogeneity amongst the constituents, a phenomenon named 
diversity-induced resonance |5j. Also, a certain degree of structural disorder in a network can 
make the propagation of information considerably more efficient [6J. 

It was recently shown that, in an ensemble of coupled excitable rotators, noise and disorder 
are able to trigger global firing During these events, a substantial fraction of the ensemble 
is excited almost synchronously, so that the joint firing of all those elements gives rise to a 
collective process which can be detected at macroscopic level. This phenomenon has been 
observed under the action of additive noise, when stochastic terms are added to the evolution 
equations, or under ensemble heterogeneity, when the parameters that define the individual 
dynamics of the rotators differ all over the ensemble. Since these previous studies considered 
all-to-all (global) homogeneous coupling between elements, the effect of disorder in their 
interaction network was not analysed. 
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Global coupling has frequently been invoked as an analytically tractable interaction pat- 
tern for which several systems are completely solvable [8j. In real situations, however, the 
interaction pattern is usually heterogeneous, with different weights for each pair of inter- 
acting elements. Some systems, such as neural networks, even require the consideration of 
interactions of different signs, to account for excitatory and inhibitory coupling [9l fTOj ITT]. 
Interaction patterns with attractive and repulsive weights in ensembles of coupled dynamical 
elements was first addressed by Daido \12 \ [T3 j fT4]. Their effect on the stability of synchro- 
nisation in an ensemble of phase oscillators, in terms of the relative number and strength 
of positive and negative couplings, has been studied recently [T5]. It was found that a cer- 
tain fraction of sufficiently strong repulsive interactions, distributed at random over the 
interaction pattern, is able to induce a transition from the highly coherent state of full 
synchronisation, where all elements oscillate in phase, to an unsynchronised state. 

The presence of randomly scattered negative coupfings in the interaction pattern of an 
ensemble of dynamical elements is a source of structural disorder, complementary to noise 
and to the heterogeneity of the ensemble. It is therefore natural to study its effects on the 
collective dynamics of coupled excitable elements and, in particular, on the phenomenon of 
global firing. This is our aim in this paper, which is organised as follows: In section [51 we 
introduce the model of coupled excitable rotators, and the order parameters used in the 
characterisation of the collective macroscopic behaviour of the ensemble. Next, in section [Sj 
we consider the case when repulsive interactions form a fully random (Erdos-Renyi) network. 
In this case, global firing is observed for finite-size systems but, due to self-averaging in the 
degree distribution of fully random networks, the phenomenon is suppressed in the limit 
of infinitely large ensembles. Thus, in section lU we introduce a model of partially random 
networks where the relative dispersion of the degree distribution does not depend on the 
size. When repulsive interactions are distributed over such a network, global firing is found 
to persist for large systems. Our conclusions are discussed in the final section [H 

2 Model and order parameters 

We consider an ensemble of coupled active rotators [16] with individual phases (j)j{t) G [0, 2tt), 
j — 1, . . . ,N, whose dynamics is given by 



The coupling strength is measured by the parameter C > 0, and each factor Wkj weights the 
interaction of a specific rotator pair. These weights are symmetric, Wkj — Wjk- Attractive 
and repulsive interactions are characterised, respectively, by Wkj > and Wkj < 0. 

The natural frequency w > is the same for all rotators. For uj > 1 and in the absence 
of coupling, the individual dynamics is oscillatory, with an actual frequency uj' = \/ uS^ — 1. 
The case lo < \ corresponds to excitable individual dynamics: for C = 0, the phase ^j of 
each rotator has two fixed points, one of them stable ((/)s < 7r/2) and the other unstable 
> 7r/2), at the two solutions of sin(/)j = w. A perturbation of the stationary stable 




(1) 



fc=i 
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solution which overcomes the unstable fixed point gives rise to the firing of the rotator, 
in the form of a long phase excursion which finally returns to the rest state (pg ■ Throughout 
this paper we focus the attention on the excitable regime u < 1. 

In our model noise is absent and, for C = 0, individual rotators are identical. Diversity 
is thus restricted to disorder in the interaction network, through the weights Wkj ■ Our main 
aim is, in fact, to analyse the effects of this source of diversity in the collective dynamics of 
the system. 

If all the phases coincide with one of the two fixed points of the individual dynamics, ei- 
ther 4>s or 4>u , the ensemble is in a stationary state -which we refer to as full synchronisation- 
for any value of C. If the weights Wkj are positive for all k and j, the fully synchronised 
state where all the rotators are in the rest state, (j)j = i/is for all j, is stable. The stability of 
this collective rest state can break down, however, in the presence of repulsive interactions, 
when some of the weights Wkj are negative and their absolute value is large enough. In this 
situation, generally, the ensemble does not reach an asymptotic stationary state. Individ- 
ual phases can now rotate irregularly around their whole domain and, as we show below, a 
regime of global firing -where a fraction of the ensemble is collectively entrained into long 
excursions from the unstable fixed point to the rest state- becomes possible. 

The collective behaviour of the rotator ensemble, including possible transitions between 
different dynamical regimes, is well characterised by a set of order parameters defined in 
terms of the individual phases 4>j(t). First, we take the average of the imaginary phase 
exponentials 

1 ^ 

p(i)exp[iM'(0] = -^exp[i(/.,(t)], (2) 

and compute the Kuramoto order parameter as p = where (•) stands for the time 

average over a long interval [5]. This parameter is a direct measure of the degree of synchro- 
nisation attained by the ensemble. For a fully synchronised state we have p — 1, whereas 
p ^ 1/Vn for a state where phases are uniformly distributed over [0, 27r). 

The Kuramoto order parameter cannot discern between the case where phases are syn- 
chronised at a fixed point, as at the rest state (f)s^ and the case where they rotate coherently, 
as expected to occur in a regime of global firing. To discriminate between static and dynamic 
entrainment, we apply the order parameter introduced by Shinomoto and Kuramoto |16j : 

C = {\p{t) exp[zvl/(f)] - {p{t) exp[zvl/(i)]) I) . (3) 

This parameter differs from zero for synchronous firing only. 

A third order parameter, frequently used in the analysis of stochastic transport [17], is 
the current 

1 ^ 

(4) 

It measures the level of (not necessarily synchronised) firing over the whole ensemble. 
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3 Erdos-Renyi networks of repulsive links 

First, we consider an ensemble of rotators governed by Eq. ([T]), whose interaction weights 
are distributed at random according to the following prescription: 

^ r 1 with probability 1 — p^, 
\ — '« with probability pd- 

The coefficient k > measures the relative strength of repulsive and attractive interactions, 
and the probability pd fixes the expected fraction of negative weights. Repulsive interactions 
define a network which, with the prescription of Eq. |[5|), has the same structure as an Erdos- 
Renyi random network [THl[l9l[2^. On the average, the network of repulsive interactions has 
PdN{N — l)/2 links. We call J\fj the set of neighbours of rotator j in that network. The mean 
number of rotators in J\fj is Pd{N — 1), and the expected dispersion around this average is 
of order V^, so that the relative dispersion decreases with the size as 1/\/N. 

It is useful to rewrite Eq. (HJ in terms of the quantities p{t) and ^'(t) introduced in Eq. 
as 

(f)j = Lu — sin (f)j + Cp sin(5' — (pj) 

+ ^ sin(0,~<^,). (6) 



N 



This equation describes the evolution of 4'j{t) as governed by the interaction with a single 
(mean-field) oscillator with phase ^{t) with an effective coupHng strength Cp{t) plus a 
negative contribution with coupling C(l -I- k) from the neighbourhood of rotator j in the 
network of repulsive links. 

3.1 Stability of the fully synchronised state 

As advanced above, in the presence of repulsive interactions, the state of full synchronisation 
may be unstable. The stability condition can be obtained from linearisation of Eq. |(6]). 
Writing 0j [t) = <j)s + 5j (t) , we get for the deviations from the fixed point the equations 

6, = -VT^s, - cs, - (-5^ - ^i)- (7) 

Full synchronisation is stable if all the eigenvalues of the matrix 

C{k + 1) 



J^-(Vl-c^2^g)I_ ^^^ ^ (M^N) (8) 

are negative or have negative real parts. Here, I is the identity matrix, M = {M^j} is the 
adjacency matrix of the network of repulsive links (defined as M^j = 1 if there is a repulsive 
link between rotators k and j, and Mkj = otherwise), and N = {Nkj} is a diagonal matrix, 
where Njj = kj = ^ ■ My is the number of rotators in Nj -namely, the degree of site j. 

Whether the stability condition for full synchronisation is fulfilled or not depends on the 
specific realisation of the network of repulsive links. Even for the same probability pd, and 
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Figure 1: Uppermost panel: Fraction fd of realisations of the network of repulsive links for which 
the fully synchronised state is unstable, as a function of the probability pd, over series of 10^ 
realisations. From left to right: = 50 (o), 100 (a) and 200 (□). The other parameters are 
C = 4, «; = 10, and u = 0.98. The vertical dotted line stands for the theoretical transition point, 
~ 9.54 X 10~^, for an infinitely large ensemble. The three lower panels show the Kuramoto 
order parameter p, the Shinomoto-Kuramoto order parameter and the current J, for the same 
realisations. 

for fixed values of C, lu and k, the stability of full synchronisation depends on the particular 
distribution of negative weights. The uppermost panel of Fig. [T] shows numerical results for 
the fraction fd of realisations of the network of repulsive links for which full synchronisation 
is unstable, as a function of pd- As expected, fd grows as the number of repulsive links 
increases. Different curves correspond to different system sizes N. 

Numerical results suggest that, for large N, there is a sudden transition from = to 
1 at a given value pjj of the probability pd- An analytical estimate of can be obtained 
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Figure 2: Left column: Time evolution of the collective phase ^'(t) (thick lines) and of the phases 
of ten representative rotators (thin lines) in an ensemble of size N = 50, for three values of the 
fraction of repulsive links: (a) pd = 0.015, (b) 0.03, and (c) 0.1. The other parameters are C = 4, 
K = 10, and w = 0.98. Right column: Time evolution of p{t), for the same realisations. The 
horizontal dashed line stands for the threshold of oscillatory behaviour, p = uj |7]. The middle 
panels correspond to the regime of global firing. 



from the assumption that the number of neighbours in the network of repulsive links is 
the same for all rotators, i.e. that the number of rotators in Afj is the same for all j. This 
approximation improves as N grows, because -as stated above- the relative dispersion in the 
number of rotators in A/} decreases as 1/ ^/N. It implies that the matrix JJ can be written as 



= -iVl-u;-' + C-CiK + l)p,)l - ^^"^^^ l 



(9) 



For iV ^ oo, the coefficient in front of the matrix M tends to zero, suggesting that its 
contribution to J can be neglected in that limit. This intuitive argument implies that the 
point at which the maximum eigenvalue becomes positive for ^ oo is 

. C 



C(l 



(10) 



A more rigorous argument, confirming this critical value for p^, can be obtained from exact 
bounds for the eigenvalues of the matrix JJ in Eq. ^ coming from the so-called semicircle 
law fST], applied to matrices M and N [l5j. The vertical dotted lines in the uppermost panel 
of Fig. [T] stands for this critical value pj. 
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Close inspection of the destabilisation of full synchronisation reveals that this transition 
is triggered by the behaviour of the rotator j* with the largest number of repulsive Hnks 
-as also observed to occur in ensembles of coupled phase oscillators [15]. For a given value 
of K, as the fraction pd of repulsive Hnks grows, the fixed point (j)s first ceases to be a 
stable state for (fij* . An arbitrarily small deviation in leads to a different equilibrium 
point, and full synchronisation breaks down, even when other rotators may remain mutually 
synchronised. For small sizes N, the value ofpd at which the rotator j* attains the number of 
negative interaction weights which makes full synchronisation unstable depends sensibly on 
the detailed structure of the network of repulsive links. This explains the smooth transition 
in fd as pd grows. For large N, on the other hand, all sites in the network of repulsive links are 
statistically equivalent with respect to the distribution of negative links. Consequently, the 
transition takes always place at the same value of pd, irrespectively of the specific realisation 
of the network. Accordingly, fd exhibits a sharp growth from to 1, which should become 
discontinuous in the thermodynamic limit. This semi-qualitative analysis sheds light on the 
role of the homogeneity of the distribution of negative Hnks in defining the nature of the 
destabilisation transition of full synchronisation. This effect of the network structure on the 
collective dynamics of the system is further studied in the foHowing. 

3.2 Numerical study of the unsynchronised regime 

In the fuHy synchronised state, the phases of aH rotators are at the stable fixed point 0s- 
The Kuramoto order parameter reaches its maximum value p = 1, whereas the absence of 
any kind of coHective motion implies that both the Shinomoto-Kuramoto parameter ( and 
the current J vanish. As full synchronisation breaks down, the parameter p is expected to 
decrease, and C and J can adopt non-zero values, as was shown in [Tj. In this section, we 
present numerical results illustrating the behaviour of the order parameters as a function of 
the fraction of repulsive links pd, for ensembles of different sizes, with C = 4, k = 10, and 
Lu = 0.98. 

Figure [2] illustrates the dynamics of the rotator ensemble for three values of the fraction 
of Pd- Results correspond to numerical integration of Eq. ^ for a system of = 50 rotators. 
In the left column, we show the time evolution of the collective phase ^{t) together with that 
of ten representative rotators. In panel (b), which corresponds to an intermediate value of 
Pd, the rotators pulse in synchrony and, consequently, the collective phase performs periodic 
rotations. We identify this behaviour with the regime of global firing. In the right column, we 
show the evolution of the parameter p{t), defined in Eq. ([2]), for the same realisations. The 
degree of synchronisation, measured by the time average of p(t), decreases as the fraction of 
repulsive interactions grows. 

In the three lower panels of Fig. [T] we plot the order parameters p, C, and J as functions of 
the fraction of repulsive links. As expected, the decay in p is initially accompanied by growth 
in both C and J. As pd increases, the current keeps growing indicating that, on average, the 
firing frequency of individual rotators is also increasing. The Shinomoto-Kuramoto parameter 
( attains a maximum at an intermediate value of and then decreases. This is an indication 
that global synchronous firing occurs when full synchronisation is unstable and for moderate 
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values oi Pd- Larger fractions of negative links, however, make the dynamical coherence of 
the population decrease, and global firing becomes less distinct. 

These results -in particular, the existence of and intermediate range of the probability 
Pd where global firing is most conspicuous- are in qualitative agreement with the collective 
behaviour of rotator ensembles where the disorder associated with repulsive interactions 
is replaced by noise and/or ensemble heterogeneity [7]. However, a noticeable quantitative 
difference resides in that, in the present case, the range of global firing becomes narrower 
as the system grows in size. The results of Fig. [T] suggest that global firing may disappear 
for sufficiently large ensembles. An explanation for this collapse can be put forward recalling 
that the instability transition of full synchronisation becomes sharper as N grows (cf. the 
uppermost panel of Fig. [TJ. As discussed at the end of section ISTTl for large systems, the 
number of repulsive links is statistically the same for all rotators. At the transition, all 
rotators abandon simultaneously the stable fixed point (ps and, consequently, the ensemble 
adopts a highly disorganised phase configuration. In small systems, on the other hand, the 
transition is gradual. The stable fixed point is first abandoned by those rotators with a large 
number of repulsive links, which can be thus entrained into global synchronous firing. 

This explanation indicates that, if it were possible to design a random network where 
the heterogeneity in the distribution of repulsive links is maintained as the system grows in 
size, the regime of global firing would persist for arbitrary large systems. In section [H we 
introduce a procedure to construct such a heterogeneous network, and numerically verify 
that the range where global firing occurs does not collapse as N grows. We also propose an 
interpolation between that kind of heterogeneous networks and Erdos-Renyi networks, and 
analyse global firing on these intermediate structures. 

4 Heterogeneous networks 

In this section, we introduce a method to construct a class of networks where the distribution 
of links per site preserves its heterogeneity as the size N increases. Specifically, the dispersion 
in the number of links per site is, in these networks, proportional to N so that the relative 
width of the distribution does not change with the size. According to the discussion at the 
end of section [3l the effects of heterogeneity observed in small random networks of repulsive 
Hnks, which disappear because of self- averaging in the random distribution of Hnks, should 
persist in these new structures. In particular, we expect to find global firing even in large 
ensembles. 

4.1 Construction of a heterogeneous network 

To build up our heterogeneous network of repulsive links, we start by taking a random 
network of size iVo and wiring probability pd- The average number of links per site is ko ~ 
PdNo, while its mean square dispersion is ctq ^ \/pdNo. Let A/^ be the set of the kj neighbours 
of site j. We now consider M identical replicas of the starting network, which we use to build 
up a network of size N = MNg. In each one of these replicas, say in replica m, we identify 
the site jm which is homologous to the site j of the starting network. The M replicas are 
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Figure 3: A heterogeneous network constructed starting by a network of A'^o = 4 sites, and 
consisting of two replicas of the starting network. Bold lines stand for links inside each replica, 
and thin lines correspond to links connecting different replicas. 

linked by connecting each site jm with all the neighbours of all its homologous sites in all 
the replicas. In other words, the site jm becomes linked to the kjM members of the sets Afj^, 
A/^M- Figure [3] illustrates the resulting network in a simple case. By construction, 
the number of neighbours of each site in the final structure equals M times the number 
neighbours of the homologous site in the starting network. Therefore, the average number of 
hnks per site and the mean square dispersion are, respectively, 

~k = MpdNo = Mfco = PdN, (11) 

and 

a = Mao = Nctq/Nq. (12) 

As advanced above, for fixed Nq, both the mean value and the dispersion in the distribution 
of links of the resulting network grow linearly with the number of sites N. The relative 
dispersion is thus independent of the network size. 

Let Mq and M be the adjacency matrices of the starting network and the resulting 
network, respectively. It is possible to show that, if Aq is an eigenvalue Mq, then A — MAq = 
NXq/Nq is an eigenvalue of M, with multiplicity Nq. The remaining eigenvalues are equal 
to zero. A similar relation holds for the matrix N of Eq. ((HI). This result differs from the 
prediction of the semicircle law, which stands for purely random matrices and states that A ~ 
VN [21]. Since, in Eq. (JH]), both M and N are divided by the network size N, the maximum 
eigenvalue of JJ -which is in turn determined by the maximum (positive) eigenvalue of M- 
does not depend on N. Thus, for fixed Nq, the critical point at which full synchronisation 
becomes unstable is independent of the size. 
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Figure 4: The order parameters as functions of the fraction of repulsive links pa, in the absence 
of rewiring, pf = 0. Different curves, which are however practically coincident, correspond to 
sizes N = 100, 200, 400, and 1000. In all cases, the network of negative links is constructed from 
a starting network of Nq = 20 sites. The other parameters are C = 4, k = 2, and uj = 0.98. 

Clearly, the network resulting from our construction is not fully random. In fact, because 
of the interconnection of all the homologous sites, it can be shown to have a large cluster- 
ing coefficient [20]. To interpolate between this structure and a fully random network, we 
introduce a mechanism of Hnk rewiring, in the spirit of small- world networks [22j • Starting 
from the heterogeneous network constructed above, we visit each site and, with probability 
Pf, rewire each of its links to a randomly chosen site all over the system. For = 1, a fully 
random network is obtained. In the following, we study the collective behaviour of an ensem- 
ble of globally coupled rotators whose network of repulsive interactions is build as described 
here, in terms of the parameters pd and pf. 

4.2 Numerical results 

In our numerical analysis, the networks of repulsive links introduced in section |4?T] are built 
starting from random networks of Nq = 20 sites, and their final sizes are N = 100, 200, 400, 
and 1000. Also, we fix C = 4, k = 2, and ui = 0.98. 
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Figure 5: The order parameters as functions of the rewiring probabihty pf, for a fixed fraction 
of repulsive hnks, pd = 0.2. Different curves correspond to sizes = 100 (o), 200 (a), 400 (o), 
and 1000 (□). In all cases, the network of negative links is constructed from a starting network 
of A^o = 20 sites. The other parameters are C = 4, k = 2, and u = 0.98. 

Figured] shows results for the order parameters p, (, and J as functions of the fraction of 
repulsive links pd, in the absence of rewiring, pf = 0. Different curves correspond to different 
sizes N. We see that, as expected, the order parameters are essentially independent of the 
system size. In particular, the regime of global firing -signalled, for intermediate values of pd, 
by the maximum in the Shinomoto-Kuramoto parameter- clearly persists for large ensembles. 

Figure [H shows the order parameters as functions of the rewiring probability pj, for a 
fixed fraction of negative links, pd = 0.2 which corresponds, approximately, to the maximum 
common firing activity in the case pj = 0, see figure HI The Shinomoto-Kuramoto parameter 
C is practically independent of for small and large p/, but exhibits a substantial dependence 
on the size for intermediate wiring probabilities. As may be expected from our results on 
Erdos-Renyi networks, in this intermediate zone, ( is larger for smaller ensembles. 

Finally, Fig. [6] summarises the dependence of the order parameters on the probabilities pd 
and Pf for an ensemble of TV = 400 rotators. The graphs represent numerical results averaged 
over series of 100 realisation for each parameter set. 
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Figure 6: The order parameters for an ensemble of = 400 rotators (A'o = 20) as functions 
of the fraction of negative hnks pd and the rewiring probabihty pf. The other parameters are 
C = 4, K = 2, and u = 0.98. 

5 Conclusion 

In this paper, we have shown that -Hke noise and ensemble heterogeneity [7]- structural 
disorder in a network of interacting excitable rotators can induce global firing of the en- 
semble. This form of collective dynamics is here a consequence of the presence of a fraction 
of repulsive couplings, randomly scattered all over the interaction network. These repulsive 
interactions destabilise the state of full synchronisation when their number and intensity are 
large enough. Just beyond this desynchronisation transition, for intermediate values of the 
fraction of repulsive couplings, the ensemble is entrained in a regime where many rotators 
are coherently excited, and global firing is triggered. While the results of this paper support 
the general mechanism for global firing put forward in [7], it is important to stress that 
the particular mechanism at work in the system studied here is quite different. In [7] the 
rotators showed heterogeneity in the form of dispersion in the natural frequencies Ui . As the 
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dispersion increased a fraction of the rotators had natural frequencies in the oscillatory range 
uJi > 1. Those rotators would spontaneously fire if they were isolated. When coupled, they 
pull the other rotators into an state of collective firing. In the system studied in this paper, 
on the other hand, the mechanism of destabilisation is produced by the effect that negative 
couplings have on individual rotators, but these would never fire if isolated. In both cases, 
however, it is the diversity, either in the natural frequencies or in the number of negative 
links, that produces the global effect. In fact, when the pattern of repulsive couplings forms 
an Erdos-Renyi random network the regime of global firing is progressively suppressed as 
the size of the ensemble grows. This effect can be ascribed to that, as the network becomes 
larger, there is self-averaging in its degree distribution, so that the number of neighbours per 
site becomes increasingly homogeneous. This homogeneity sharpens the desynchronisation 
transition and, at the same time, causes the collapse of the zone where global firing is pos- 
sible for smaller systems. To avoid this effect of self-averaging -which clearly illustrates how 
increasing homogeneity can inhibit certain forms of coherent behaviour- we have introduced 
an algorithm to construct partially random networks whose degree distribution maintains its 
relative dispersion as the network grows. In this case, in fact, the order parameters which 
characterise the macroscopic behaviour of the rotator ensemble become essentially indepen- 
dent of the size. In particular, the regime of global firing persists for large systems. 

The present results extend previous work on the effects of diversity of various origins over 
the collective dynamics of complex systems. These forms of diversity include now heteroge- 
neous interaction patterns, with positive (attractive) and negative (repulsive) coupHngs. The 
extension is thus relevant to such systems as neural networks, where interactions of different 
signs are present in the form of activator and inhibitory synapses [lOj. 
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